Accurate determination of crystal structures based on averaged local bond order 

parameters 
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Local bond order parameters based on spherical harmonics, also known as Steinhardt order pa- 
rameters, are often used to determine crystal structures in molecular simulations. Here we propose 
a modification of this method in which the complex bond order vectors are averaged over the first 
neighbor shell of a given particle and the particle itself. As demonstrated using soft particle systems, 
this averaging procedure considerably improves the accuracy with which different crystal structures 
can be distinguished. 



I. INTRODUCTION 

In computational studies of crystallization from the 
undercooled liquid one needs to be able to distinguish 
particles that are part of the crystal from those that be- 
long to the liquid. Ideally, such an assignment should be 
based on the local environment of the particles only. One 
method to do that, which is independent of the specific 
crystal structure and does not require the definition of a 
reference frame, is provided by the following algorithm 
based on spherical harmonics First, the complex vec- 
tor qi m (i) of particle i is defined as 



qim(i) 



1 



N b (i) 



N b (i) 
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Here, Nb(i) is the number of nearest neighbors of particle 
i, I is a free integer parameter and m is an integer that 
runs from m = — I to m = +?. The functions Yi m {jCij) 
arc the spherical harmonics and is the vector from 
particle i to particle j. Using the set of complex vectors 
<?6m one can then define the scalar product 



(2) 



which measures the correlation between the structures 
surrounding particles i and j. The * indicates complex 
conjugation. Two particles i and j are defined to be 
connected if Sij is larger than a given value, typically 
Sij > 0.5. A particle is solid-like if the number of connec- 
tions it has with its neighbors is above a certain thresh- 
old, typically between 6 and 8. If a particle is connected 
to less particles, it is considered to be liquid-like. Using 
this criterion to distinguish solid-like from liquid- like par- 
ticles one can then search for clusters of connected solid- 
like particles. The size of the largest of these crystalline 
clusters N is often used as a reaction coordinate to follow 
the progress of the crystallization process. Provided that 
this reaction coordinate captures the essential physics of 
the crystal nucleation, the Gibbs free energy G(N) cal- 
culated as a function of N provides a means to estimate 
the rate at which the nucleation of the crystalline phase 
occurs. 



This procedure very efficiently distinguishes between 
solid-like and liquid-like particles, but does not discrim- 
inate between different crystal structures. A set of pa- 
rameters which hold the information of the local struc- 
ture are the local bond order parameters, or Steinhardt 
order parameters, defined as [2| 
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Depending on the choice of I, these parameters are sen- 
sitive to different crystal symmetries. Each of them de- 
pends on the angles between the vectors to the neigh- 
boring particles only and therefore these parameters are 
independent of a reference frame. Different approaches 
based on these local bond order parameters were devel- 
oped to analyze the structure of the crystalline nucleus 
during the freezing event. Especially q± and qe are often 
used as they are a good choice to distinguish between 
cubic and hexagonal structures 0, 0, HI ■ 

In practice, local bond order parameters are used in 
different ways. Frenkcl and coworkers 0, @, 0, 3 ana- 
lyzed the structure of crystalline clusters in terms of or- 
der parameter distributions. To do that, they first com- 
puted the distributions of q^ and q§ for the liquid and for 
perfect FCC and BCC crystals. Due to thermal fluctua- 
tions, these distributions can be rather broad. Then, they 
determined the distributions of the same order parame- 
ters in the crystalline cluster. These distributions are 
represented as a superposition of the distribution func- 
tions of the perfect phases. The superposition coefficients 
cfcc, cbcc and cliq, determinded by mean square min- 
imization, then yield information on the composition of 
the cluster. For instance, cfcc ~ 0.5, cbcc ~ 0.5, and 
c liq ~ would be indicative of a cluster that is half in 
the FCC structure and half in the BCC structure. 

Another bond order parameter method, used for in- 
stance in Rcf. [9f], defines regions in the two dimensional 
g4-gg-plane and g4-W4-planc. The crystalline structure 
around a given particle is characterized by its positions 
in these planes. The parameter wi necessary for that 
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analysis is defined as 



parameter wi, 
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Here, the integers mi, 171% and m3 run from — I to +1, but 
only combinations with mi + vcii + = are allowed. 
The term in brackets is the Wigncr 3-j symbol [lfj • Using 
this approach, one can determine the type of crystalline 
structure occurring around each individual particle. 

As mentioned above, thermal fluctuations smear out 
the order parameter distributions such that it may be 
difficult to distinguish local crystalline structures based 
on Steinhardt bond order parameters. In the following 
we present a simple method to increase the accuracy 
of the crystal structure determination by averaging over 
the bond order parameters of nearest neighbor particles. 
This averaging procedure is discussed in Sec. |n]and val- 
idated for two different soft sphere systems in Sec. IIIII 
Some conclusions are provided in Sec. IIVI 



II. AVERAGED BOND ORDER PARAMETERS 

The crystal structure determination described above 
can be improved by using the following averaged form of 
the local bond order parameters: 



where 
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Here, the sum from k = to Nb(i) runs over all neighbors 
of particle i plus the particle i itself. Thus, to calculate 
q~l(i) of particle i one uses the local orientational order 
vectors qi m (j-) averaged over particle i and its surround- 
ings. While qi(i) holds the information of the structure 
of the first shell around particle i, its averaged version 
q~i(i) also takes into account the second shell. One might 
say that using the parameter q\ instead of qi increases 
the accuracy of the distinction of different structures at 
the price of a coarsening of the spacial resolution. In this 
sense, the averaged local bond order parameter is similar 
to the scalar product of Eq. © used to decide whether a 
particle is in a solid- like or liquid- like environment. Also 
in that case the second particle shell is effectively taken 
into account. 

The averaged orientational order parameter qi m can also 
be used to define an averaged version wi of the order 
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III. RESULTS 



(7) 



In the following, we verify that these averaged forms 
of the bond order parameters indeed increase the accu- 
racy of the crystal structure determination. The calcu- 
lations are done for two distinct systems: one in which 
the particles interact via a Lennard- Jones potential and 
one in which the interaction is of the Gaussian core form 

Pimm. 



A. Lennard-Jones system 

For the calculations of the Lennard-Jones system [l4| 
the temperature was k^T = 0.92e and the pressure was 
P = 5. 68eer -3 , corresponding to 20% undercooling of the 
liquid phase [6(. The same conditions were used in Ref. 
[J to study homogeneous crystal nucleation. This phase 
point corresponds to a mean density of p ~ 1.05<r in 
the FCC crystal, which is the most stable phase un- 
der these conditions. All simulations were carried out 
in the isobaric-isothermal ensemble, but calculations at 
constant volume for the corresponding densities yielded 
essentially identical results. The particle number was 
N = 432 'in the BCC crystal, N = 864 in the FCC crys- 
tal, and N = 512 in the HCP crystal and in the under- 
cooled liquid. Two particles were defined to be neighbors 
if their distance was smaller than r/v = \Ao. 

In Fig. Q] we compare the distributions of the local 
bond order parameters q<± and q§ with those of the aver- 
aged bond order parameters q~4 and qg. Due to the av- 
eraging procedure, the overlap between the distributions 
decreases considerably thus allowing one to distinguish 
between different structures more precisely. As can be 
inferred from Fig. [1] the reduced overlap is due not only 
to a narrowing of the distributions but also to a shift 
of the means towards smaller values. While one would 
expect the order parameter distributions to narrow on 
the basis of the central limit theorem, one is tempted 
to attribute the shift of the means to local correlations. 
To test this speculation, we have artificially eliminated 
correlations by computing qi from uncorrelated values of 
qi m generated with the appropriate statistics. Also in 
this case the narrowing of the distributions was accom- 
panied by a shift of the means indicating that this shift 
is not due to local correlations between the local order 
parameters of nearest neighbors but rather originates in 
the functional form of qi . 
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FIG. 1: Top: Probability distributions of 94 (solid lines) and 
q4 (dashed lines) for the FCC, BCC and HCP crystals and for 
the undercooled liquid (LIQ) in the Lennard- Jones system. 
Bottom: Probability distributions of q& (solid lines) and q§ 
(dashed lines) for the same phases. 





94 <?4 


96 96 


BCC 
FCC 
HCP 
LIQ 


0.089988 0.033406 
0.170880 0.158180 
0.107923 0.084052 
0.109049 0.031246 


0.440526 0.408018 
0.507298 0.491385 
0.445384 0.421813 
0.360012 0.161962 



TABLE I: Mean of the distributions of g4, 94, ge, and 96 for 
the BCC, FCC, and HCP crystals and the undercooled liquid 
in the Lennard- Jones system. 

Tableland Tablc[Il]show the mean value and the stan- 
dard deviation of the distributions of (74, (74, qe and qe, 
respectively, for all four phases we studied. The stan- 
dard deviation decreases by a factor of 2-4 in most cases. 
To quantify the improvement of the averaged order pa- 
rameter we define the overlap between two distributions 



as 

O a p = -!- / P a {x)P p {x)dx. (8) 

-<v Q( 3 J 

Here, the indices a and (3 denote the various crystal struc- 
tures and N a p is given by 



N a0 = Jj PZ(x)dx J P*(y)dy. (9) 





g4 g4 


g6 g6 


BCC-FCC 
BCC-HCP 
BCC-LIQ 
FCC-HCP 
FCC-LIQ 
HCP-LIQ 


0.124377 0.000000 
0.738776 0.007880 
0.881215 0.992978 
0.210924 0.000188 
0.233250 0.000000 
0.938029 0.001096 


0.33959496 0.017147 
0.94732246 0.866026 
0.46374455 0.000004 
0.31311693 0.033288 
0.14914436 0.000000 
0.35744697 0.000000 



TABLE III: Overlap between the distributions of g4, g4, g6, 
and g6 for the various phases in the Lennard- Jones system. 

In Table irni we summarize the overlap between the dis- 
tributions of (74, 94, qe, and q@ for BCC, FCC, and HCP 
crystals and the undercooled liquid. In most cases the 
overlap decreases by several orders of magnitude when 
going from the local order parameter to the averaged ver- 
sion. One exception are the (74 distributions of the liquid 
and the BCC crystal, which have an overlap of nearly 1 
even for the averaged order parameter. All phases are, 
however, resolved very well using two averaged bond or- 
der parameters. Figure [2] shows a comparison of the var- 
ious crystalline and liquid phases as scatter plots in the 
<74-g6-pla n e, the q4-qe-plsLne, the <74-u>4-plane, and the W4- 
u>4-plane. For the averaged version of the order param- 
eters (right) the four phases separate much better. In 
practice, the separation between BCC and HCP is the 
most difficult one. In the qi-qe-plane these two crystal 
structures can be distinguished easily and also in the yV 
u>4-plane they are well separated. When going from q 4 
to 54 the overlap between these two structures decreases 
by two orders of magnitude (see Tab. IIII[) . 

B. Gaussian Core Model 





g4 g4 


g6 g6 


BCC 
FCC 
HCP 
LIQ 


0.026831 0.010782 
0.032787 0.014346 
0.019476 0.009434 
0.031992 0.008786 


0.034791 0.020516 
0.043301 0.020566 
0.028992 0.015965 
0.066518 0.039360 



TABLE II: Standard deviation of the distributions of g4, 54, 
g6, and g6 for the BCC, FCC, and HCP crystals and the 
undercooled liquid in the Lennard- Jones system. 



For the Gaussian core system, in which particles in- 
teract via the pair-interaction v(r) = e exp[— (r/cr) 2 ] 
[ill . H2I HU, we chose a temperature of fc^T = 0.0033e 
and pressure P — 0.01 leer— 3 at which the thermody- 
namically most stable phase is an FCC crystal with den- 
sity O.llcr" 3 . These conditions correspond to the same 
level of undercooling as that of the Lennard-Jones sys- 
tem described above. All calculations were carried out in 
the isothermal-isobaric ensemble and the particle num- 
bers for the various crystal structures where the same 
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FIG. 2: Top: Comparison between the g4-96-plane (left) and 
the (?4-q6-plane (right) for the Lennard- Jones system in three 
different crystalline structures and in the liquid phase. Each 
point corresponds to a particular particle, where 2000 points 
from each structure were chosen randomly. Bottom: 54-^4- 
plane (left) and the g4-«)4-plane (right). 



as in the Lennard-Jones system. The nearest neighbor 
distance was rpf = 3.0a. 

Also for the Gaussian core system we find a shift in 
the mean and a decrease in the width of the distributions 
(see Fig. [3J . The mean and the standard deviation of the 
distributions for the various phases are listed in Tables 
HVI andlVl respectively. The overlap between the different 
structures is shown in Table IVIl In most cases the overlap 
decreases by a large factor, only the overlap in q± between 
BCC and liquid increases. But as in the Lennard-Jones 
case all structures can be distinguished very well if two 
order parameters are taken into account (see Fig. 





94 94 


96 96 


BCC 
FCC 
HCP 
LIQ 


0.085581 0.031728 
0.155336 0.134388 
0.109723 0.073369 
0.126950 0.040297 


0.437129 0.407515 
0.474079 0.447782 
0.424627 0.385720 
0.375121 0.158913 



TABLE IV: Mean of the distributions of 94, 54, 96, and 95 for 
the BCC, FCC, and HCP crystals and the undercooled liquid 
in the Gaussian core system. 



IV. CONCLUSIONS 

The averaged bond order parameters proposed in this 
paper separate different crystal structures more accu- 
rately than regular bond order parameters. Instead of 
using the first shell of surrounding particles only, this 
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FIG. 3: Top: Probability distributions of 94 (solid lines) and 
94 (dashed lines) for the FCC, BCC and HCP crystals and 
for the undercooled liquid (LIQ) in the Gaussian core system. 
Bottom: Probability distributions of 96 (solid lines) and 96 
(dashed lines) for the same phases. 





94 94 


96 96 


BCC 
FCC 
HCP 
LIQ 


0.023841 0.008024 
0.039924 0.018200 
0.023617 0.011605 
0.037891 0.011282 


0.037274 0.021006 
0.053852 0.027126 
0.038462 0.022251 
0.062874 0.038372 



TABLE V: Standard deviation a of the distributions of §4, 
94, 96, and qe for the BCC, FCC, and HCP crystals and the 
undercooled liquid in the Gaussian core system. 



method also takes into account the second shell, just as 
one usually does when distinguishing between solid-like 
and liquid- like particles (see Eq. ([2])). Due to the averag- 
ing procedure, the overlap between the order parameter 
distributions belonging to different phases decreases by 
up to a few orders of magnitude. The sharper distinction 
between the different phases is obtained at the cost of a 
slightly reduced spatial resolution. Wc demonstrated the 
enhancement in the crystal structure determination for 
two systems of soft particles interacting via a Lennard- 
Jones potential and a Gaussian core potential. For the 
same degree of undercooling the improvement is similar 
in both cases. Particularly in the q^-qQ-pl&rie the sepa- 
ration of the liquid phase from the solid phases is pro- 
nounced indicating that it might be possible to use this 
method also to distinguish between liquid and solid parti- 
cles if only BCC, FCC and HCP structures are expected 
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§4 Q4 


qe qe 


BCC-FCC 
BCC-HCP 
BCC-LIQ 
FCC-HCP 
FCC-LIQ 
HCP-LIQ 


0.248341 0.000007 
0.687131 0.015436 
0.619532 0.832443 
0.496870 0.022218 
0.632466 0.000261 
0.928871 0.128086 


0.685253 0.462343 
0.964792 0.753398 
0.581065 0.000006 
0.624299 0.197885 
0.445508 0.000000 
0.682437 0.000018 



TABLE VI: Overlap between the distributions of 94, Q4, qe, 
and <j6 for the various phases in the Gaussian core system. 



to be important for the process under study. We expect 
that the averaging procedure described in this paper lead 
to similar improvements also for other Steinhardt bond 
order parameters, e.g., q$ and w$. 




FIG. 4: Top: Comparison between the g4-Q6-plane (left) and 
the g4-Q6-plane (right) in the Gaussian core system . Bottom: 
Q4-u;4-plane (left) and the g4-W4-plane (right). 
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